!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!> \brief
!! Flow past a circular cylinder.
!!
!! \n
!!
!! The initial condition is a uniform flow; we enforce Dirichlet
!! boundary on the inflow domain and homogeneous Neumann conditions on
!! the remaining boundaries. Boundary markers are set as follows:
!! <ul>
!!   <li> 3 for the inflow region, where we have
!!   \f$\underline{u}_D=U\underline{e}_1\f$, for some
!!   \f$U\in\mathbb{R}\f$ and with \f$\underline{e}_1\f$ denoting the
!!   unit vector along the first coordinate axis;
!!   <li> values from 5 to 8 for the cylinder, where we enforce
!!   \f$\underline{u}_D=0\f$;
!!   <li> any other value for the remaining boundary regions,
!!   where we use homogeneous Neumann boundary conditions.
!! </ul>
!!
!! There is no known analytic solution for this test case; some
!! reference solutions can be found in <a
!! href="http://dx.doi.org/10.1063/1.868939">[Henderson, Barkley,
!! 1996]</a>, <a
!! href="http://journals.cambridge.org/action/displayAbstract?fromPage=online&aid=13139&fulltextType=RA&fileId=S0022112096004326">[Prasad,
!! Williamson, 1997]</a>, <a
!! href="http://dx.doi.org/10.1063/1.1338544">[Wen, Lin, 2001]</a>, <a
!! href="http://dx.doi.org/10.1002/fld.486">[Snyder, Degrez, 2003]</a>
!! and <a href="http://dx.doi.org/10.1063/1.2844875">[Inoue, Akira,
!! 2008]</a>.  As reported in <a
!! href="http://dx.doi.org/10.1063/1.1338544">[Wen, Lin, 2001]</a>, in
!! the range \f$46<Re<1000\f$ the Strouhal and Reynolds numbers are
!! related by
!! \f{displaymath}{
!!   St = 0.2417 - 0.8328\,Re^{-0.4808}e^{-0.001859\,Re}
!! \f}
!! where the two dimensionless numbers are
!! \f{displaymath}{
!!   Re = \frac{U_\infty d}{\nu}, \qquad St = \frac{fd}{U_\infty}
!! \f}
!! and the symbols indicate the asymptotic velocity \f$U_\infty\f$,
!! the cylinder diameter \f$d\f$, the kinematic viscosity \f$\nu\f$
!! and the vortex frequency \f$f\f$.
!<----------------------------------------------------------------------
module mod_cylinder_test

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_fu_manager, only: &
   mod_fu_manager_initialized, &
   new_file_unit

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_cylinder_test_constructor, &
   mod_cylinder_test_destructor,  &
   mod_cylinder_test_initialized, &
   test_name, &
   test_description,&
   coeff_visc,&
   coeff_f,   &
   coeff_dir, &
   coeff_init

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members
 character(len=*), parameter ::    &
   test_name = "cylinder"

! Module variables

 ! public members
 character(len=100), protected ::    &
   test_description(3)
 logical, protected ::               &
   mod_cylinder_test_initialized = .false.
 ! private members
 real(wp) :: nu, ub ! viscosity and imposed velocity
 character(len=*), parameter :: &
   test_input_file_name = 'cylinder_test.in'
 character(len=*), parameter :: &
   this_mod_name = 'mod_cylinder_test'

 ! Input namelist
 namelist /input/ &
   nu, ub

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_cylinder_test_constructor(d)
  integer, intent(in) :: d

  integer :: fu, ierr
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.)    .or. &
       (mod_fu_manager_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_cylinder_test_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   ! Read input file
   call new_file_unit(fu,ierr)
   open(fu,file=trim(test_input_file_name), &
      status='old',action='read',form='formatted',iostat=ierr)
    if(ierr.ne.0) call error(this_sub_name,this_mod_name, &
      'Problems opening the input file')
    read(fu,input)
   close(fu,iostat=ierr)

   write(test_description(1),'(a,i1)') &
     'Cavity test case, dimension d = ',d
   write(test_description(2),'(a,e9.3,a)') &
     '  viscosity nu = ',nu,';'
   write(test_description(3),'(a,e9.3,a)') &
     '  boundary velocity Ub = ',ub,'.'

   mod_cylinder_test_initialized = .true.
 end subroutine mod_cylinder_test_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_cylinder_test_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_cylinder_test_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_cylinder_test_initialized = .false.
 end subroutine mod_cylinder_test_destructor

!-----------------------------------------------------------------------
 
 pure function coeff_visc(x) result(nu_x)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: nu_x(size(x,2))

   nu_x = nu
 end function coeff_visc
 
!-----------------------------------------------------------------------
 
 pure function coeff_f(x) result(f)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: f(size(x,1),size(x,2))

   f = 0.0_wp
 end function coeff_f
 
!-----------------------------------------------------------------------

 pure function coeff_dir(x,breg) result(d)
  real(wp), intent(in) :: x(:,:)
  integer, intent(in) :: breg
  real(wp) :: d(size(x,1),size(x,2))

   d = 0.0_wp
   if(breg.eq.3) d(1,:) = ub
 
 end function coeff_dir
 
!-----------------------------------------------------------------------

 pure function coeff_init(x) result(u)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: u(size(x,1),size(x,2))

   ! take an initial angle to break symmetry and simplify the onset of
   ! the alternating flow
   u(1,:) = ub
   u(2,:) = tan(-30.0_wp*3.14_wp/180.0_wp)*ub
   if(size(u,1).gt.2) u(3:,:) = 0.0_wp
 
 end function coeff_init
 
!-----------------------------------------------------------------------

end module mod_cylinder_test

